Single Cell RNA-Seq: Integration & Visualization

Load packkages

UMAP annotated clusters

scRNA_epi_saline <- readRDS("/Users/alexanderfoote/dev/Projs/R/Upper_Airway_Coding_Project/rds/scRNA_epi_saline_rpca_refined_annotated_11122024.rds")

#Assign annotations for larynx vs tracheal regions
DefaultAssay(scRNA_epi_saline) <- "RNA"
Idents(scRNA_epi_saline) <- "seurat_clusters"
scRNA_epi_saline <- RenameIdents(scRNA_epi_saline, 
                   '0' =  'Tracheobronchial',
                  '1' =  'Tracheobronchial',
                  '2' =  'Tracheobronchial',
                  '3' =  'Tracheobronchial',
                  '4' =  'Tracheobronchial',
                  '11' =  'Tracheobronchial',
                  '18' = 'Tracheobronchial',
                   '15' = 'Tracheobronchial',
                  '16' =  'Submucosal-gland',
                  '9' = 'Submucosal-gland',
                  '13' = 'Submucosal-gland',
                  '7' = 'Submucosal-gland',
                  '10' = 'Submucosal-gland',
                  '5' = 'Submucosal-gland',
                  '12' =   'Pharyngolaryngeal',
                  '8' =   'Pharyngolaryngeal',
                  '14' =   'Pharyngolaryngeal',
                  '6' =   'Tracheobronchial',
                  '17' =  'Pharyngolaryngeal')

name <- "scRNA_epi_labels"
plot2 <- DimPlot(scRNA_epi_saline,reduction = "umap", pt.size=.1) + NoLegend()
LabelClusters(plot2, id="ident", size=5,repel=T, box.padding=.5)

#Assign cluster annotation for all cell types
scRNA_epi_saline$Region <- Idents(scRNA_epi_saline)

#Assign annotations to seurat_clusters
DefaultAssay(scRNA_epi_saline) <- "RNA"
Idents(scRNA_epi_saline) <- "seurat_clusters"
scRNA_epi_saline <- RenameIdents(scRNA_epi_saline, 
                  '0' =  'Basal-v.trachea',
                  '1' =  'Basal-d.trachea',
                  '12' =  'Basal-larynx',
                  '8' =  'Parabasal-larynx',
                  '14' =  'Suprabasal-larynx',
                  '6' =  'Secretory-trachea',
                  '17' =  'Secretory-larynx',
                  '16' =  'Basal-myoepithelial',
                  '2' =  'Club-proximal',
                  '3' =  'Club-mid',
                  '4' =  'Club-distal',
                  '11' = 'Club-proximal',
                  '9' =  'Serous-duct',
                  '5' = 'Goblet-2',
                  '10' = 'Serous-acini',
                  '7' = 'Goblet-2',
                  '13' = 'Goblet-1',
                  '18' =  'Tuft/NE/Ionocyte',
                  '15' = 'Ciliated')

name <- "scRNA_epi_labels"
plot1 <- DimPlot(scRNA_epi_saline,reduction = "umap", pt.size=.1) + NoLegend()
LabelClusters(plot1, id="ident", size=5,repel=T, box.padding=.5)

#Assign cluster annotation for all cell types
scRNA_epi_saline$CellType <- Idents(scRNA_epi_saline)
LabelClusters(plot1, id="ident", size=5,repel=T, box.padding=.5)

LabelClusters(plot2, id="ident", size=5,repel=T, box.padding=.5)

Add custom cluster annotation

NE <- WhichCells(object = scRNA_epi_saline, expression = Ascl1 > 1)
parentcluster <- WhichCells(object = scRNA_epi_saline, idents = 'Tuft/NE/Ionocyte')
NE <- NE[NE %in% parentcluster]
Idents(scRNA_epi_saline, cells = NE) <- 'Neuroendocrine'

Tuft <- WhichCells(object = scRNA_epi_saline, expression = Trpm5 > 1)
parentcluster <- WhichCells(object = scRNA_epi_saline, idents = 'Tuft/NE/Ionocyte')
Tuft <- Tuft[Tuft %in% parentcluster]
Idents(scRNA_epi_saline, cells = Tuft) <- 'Tuft'

#Reassign cluster annotation
scRNA_epi_saline$CellType <- Idents(scRNA_epi_saline)

scRNA_epi_saline$CellType <- factor(scRNA_epi_saline$CellType,levels=c('Basal-d.trachea','Basal-v.trachea','Basal-larynx','Parabasal-larynx','Basal-myoepithelial','Club-proximal','Club-mid','Club-distal',"Serous-duct",'Serous-acini','Goblet-1','Goblet-2','Secretory-unknown','Secretory-larynx','Suprabasal-larynx',"Neuroendocrine",'Tuft','Ciliated'))

scRNA_epi_saline$CellType <- factor(scRNA_epi_saline$CellType,levels=c('Basal-trachea','Basal-larynx','Parabasal-larynx','Basal-myoepithelial','Club-proximal','Club-mid','Club-distal',"Serous-duct",'Serous-acini','Goblet-1','Goblet-2','Secretory-trachea','Secretory-larynx','Suprabasal-larynx',"Neuroendocrine",'Tuft','Ciliated'))


#Remove 'Tuft/NE/Ionocyte' cluster
DimPlot(scRNA_epi_saline,reduction = "umap",label = TRUE,pt.size=1) 

#Reassign cluster annotation
scRNA_epi_saline$CellType <- Idents(scRNA_epi_saline)

name <- "scRNA_epi_labels_custom-annotate_celltype"
DimPlot_scCustom(scRNA_epi_saline,reduction = "umap",label = T,pt.size=.3, DiscretePalette_scCustomize(num_colors = 24, palette = "alphabet"), figure_plot = TRUE, label.size = 4, label.box = T) 

name <- "scRNA_epi_labels_custom-annotate_VlnPlotQC"
VlnPlot(scRNA_epi_saline, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3, pt.size = 0)

Visualize all marker genes on annotated seurat clusters per cell type

#Per Cell Type
DefaultAssay(scRNA_epi_saline) <- "RNA"
Idents(scRNA_epi_saline) <- "CellType"
scRNA_epi_saline$CellType <- factor(scRNA_epi_saline$CellType,levels=c('Basal-d.trachea','Basal-v.trachea','Basal-larynx','Parabasal-larynx','Basal-myoepithelial','Suprabasal-larynx','Secretory-trachea','Secretory-larynx','Club-proximal','Club-mid','Club-distal','Goblet-1','Goblet-2',"Serous-duct",'Serous-acini',"Neuroendocrine",'Tuft','Ciliated'))
celltype_markers <- c("Trp63","Cav1","Tgm2","Igfbp2","Tmprss11a","Acta2","Krt13","Kcnj16","Duoxa2","Slc6a15","Ces1f","Scgb1a1","Tff2","Lipf","Slc34a2","Bpifb5","Ascl1","Trpm5","Foxj1")
DotPlot_scCustom(scRNA_epi_saline, features=celltype_markers, colors_use= c("#E4E1E3FF","#66B0FFFF","#F6222EFF"), flip_axes = F, x_lab_rotate = TRUE)

# Find markers and limit to those expressed in greater than 60% of target population
all_markers <- FindAllMarkers(object = scRNA_epi_saline) %>%
    Add_Pct_Diff() %>%
    filter(pct_diff > 0.6)
top_markers <- Extract_Top_Markers(marker_dataframe = all_markers, num_genes = 5, named_vector = FALSE,
    make_unique = TRUE)
Clustered_DotPlot(scRNA_epi_saline, features = top_markers, k = 14)

## [[1]]

## 
## [[2]]

DefaultAssay(scRNA_epi_saline) <- "integrated"
DoHeatmap(scRNA_epi_saline, angle = 90,size = 3, group.by = "CellType", features = celltype_markers) + scale_fill_gradientn(colors = c("#E4E1E3FF","#66B0FFFF","#F6222EFF"))

#Per Cell Type
DefaultAssay(scRNA_epi_saline) <- "RNA"
Idents(scRNA_epi_saline) <- "CellType"
scRNA_epi_saline$CellType <- factor(scRNA_epi_saline$CellType,levels=c('Basal-d.trachea','Basal-v.trachea','Basal-larynx','Parabasal-larynx','Basal-myoepithelial','Suprabasal-larynx','Secretory-trachea','Secretory-larynx','Club-proximal','Club-mid','Club-distal','Goblet-1','Goblet-2','Serous-acini',"Serous-duct","Neuroendocrine",'Tuft','Ciliated'))
celltype_markers <- c("Trp63","Krt5","Krt17","Cav1","Tgm2","Dapk1","Lgr6","Igfbp2","Tmprss11a","Ntng1","Krt14","Ntrk3","Cxcl12","Acta2","Cxcl14","Krt13","Krt6a","Tmprss11b","Rbp2","Sprr1a","Muc1","Muc4","Muc20","Tnfaip2","Kcnj16","Tnfsf10","Il1a","Muc13","Bpifa1","Scgb1a1","Scgb3a2","Cyp2a5","Scgb3a1","Muc5b","Tff2","Lman1l","Wfdc12","Lipf","Dcpp3","Ltf","Dmbt1","Bpifb1","Bpifb5","Lyz1","Clu","Cryab","Slc34a2","Ascl1","Ngf","Cxcl13","Calca","Trpm5","Gnat3","Pou2f3","Lrmp","Dclk1","Spib","Foxj1","Ccdc39","Tmem212","Dynlrb2")
DotPlot_scCustom(scRNA_epi_saline, features=celltype_markers, colors_use= c("#E4E1E3FF","#66B0FFFF","#F6222EFF"), flip_axes = F, x_lab_rotate = TRUE)

DefaultAssay(scRNA_epi_saline) <- "integrated"
DoHeatmap(scRNA_epi_saline, angle = 90,size = 3, group.by = "CellType", features = celltype_markers) + scale_fill_gradientn(colors = c("#E4E1E3FF","#66B0FFFF","#F6222EFF"))

Highlight a subset of cells

# Universal Step 1: List all identities and confirm they match the intended clusters
unique_clusters <- unique(Idents(scRNA_epi_saline))  # Check if this includes the desired identities

# Step 2: Create a named color vector FOR BASAL CELLS where each identity is explicitly matched to a color
# Assign grey to unhighlighted clusters by default
colors_use_basal <- rep("grey", length(unique_clusters))
names(colors_use_basal) <- unique_clusters

# Manually set colors for BASAL highlighted cell types
colors_use_basal["Basal-larynx"] <- "#993F00FF"
colors_use_basal["Basal-v.trachea"] <- "#0075DCFF"
colors_use_basal["Basal-d.trachea"] <- "#F0A0FFFF"
colors_use_basal["Basal-myoepithelial"] <- "#191919FF"

colors_use_basal["Basal-larynx"] <- "#993F00FF"
colors_use_basal["Basal-trachea"] <- "#F0A0FFFF"
colors_use_basal["Basal-myoepithelial"] <- "#191919FF"

# Step 3: Use DimPlot with this named color vector
# Re-order cell types
scRNA_epi_saline$CellType <- factor(scRNA_epi_saline$CellType,levels=c('Basal-v.trachea','Basal-d.trachea','Basal-larynx','Parabasal-larynx','Basal-myoepithelial','Suprabasal-larynx','Secretory-trachea','Secretory-larynx','Club-proximal','Club-mid','Club-distal','Goblet-1','Goblet-2',"Serous-duct",'Serous-acini',"Neuroendocrine",'Tuft','Ciliated'))

Idents(scRNA_epi_saline) <- "CellType"
DimPlot(
  scRNA_epi_saline, 
  group.by = "CellType",   # Ensure "CellType" is the correct metadata for grouping
  reduction = "umap",
  label = F, 
  pt.size = 1, 
  cols = colors_use_basal  # Named vector applies specific colors to specific identities
)

VlnPlot(
  scRNA_epi_saline, 
  features = c("Trp63"), 
  group.by = "CellType", 
  cols = colors_use_basal  # Set custom colors
)

# Step 2: Create a named color vector FOR SMG where each identity is explicitly matched to a color
# Assign grey to unhighlighted clusters by default
colors_use_SMG <- rep("grey", length(unique_clusters))
names(colors_use_SMG) <- unique_clusters

# Manually set colors for SMG highlighted cell types
colors_use_SMG["Serous-duct"] <- "#003380FF"
colors_use_SMG["Serous-acini"] <- "#FFA405FF"
colors_use_SMG["Goblet-1"] <- "#9DCC00FF"
colors_use_SMG["Goblet-2"] <- "#C20088FF"
colors_use_SMG["Basal-myoepithelial"] <- "#191919FF"

# Step 3: Use DimPlot with this named color vector
DimPlot(
  scRNA_epi_saline, 
  group.by = "CellType",   # Ensure "CellType" is the correct metadata for grouping
  reduction = "umap",
  label = F, 
  pt.size = 1, 
  cols = colors_use_SMG  # Named vector applies specific colors to specific identities
)

# Step 2: Create a named color vector FOR TERMINAL CELLS where each identity is explicitly matched to a color
# Assign grey to unhighlighted clusters by default
colors_use_ter <- rep("grey", length(unique_clusters))
names(colors_use_ter) <- unique_clusters

# Manually set colors for SMG highlighted cell types
colors_use_ter["Club-proximal"] <- "#808080FF" 
colors_use_ter["Club-mid"] <- "#94FFB5FF"
colors_use_ter["Club-distal"] <-  "#8F7C00FF"
colors_use_ter["Ciliated"] <- "#FF0010FF" 
colors_use_ter["Neuroendocrine"] <- "#FFA8BBFF"
colors_use_ter["Tuft"] <- "#005C31FF"
colors_use_ter["Suprabasal-larynx"] <- "#426600FF"
colors_use_ter["Secretory-trachea"] <- "#2BCE48FF"
colors_use_ter["Secretory-larynx"] <- "#FFCC99FF" 

# Step 3: Use DimPlot with this named color vector
DimPlot(
  scRNA_epi_saline, 
  group.by = "CellType",   # Ensure "CellType" is the correct metadata for grouping
  reduction = "umap",
  label = F, 
  pt.size = 1, 
  cols = colors_use_ter  # Named vector applies specific colors to specific identities
)

# Step 2: Create a named color vector FOR LARYNX CELLS where each identity is explicitly matched to a color
# Assign grey to unhighlighted clusters by default
colors_use_larynx <- rep("grey", length(unique_clusters))
names(colors_use_larynx) <- unique_clusters

# Manually set colors for SMG highlighted cell types
colors_use_larynx["Basal-larynx"] <- "#993F00FF"
colors_use_larynx["Parabasal-larynx"] <- "#FFCC99FF" 
colors_use_larynx["Suprabasal-larynx"] <- "#4C005CFF" 
colors_use_larynx["Secretory-larynx"] <- "#005C31FF"

# Step 3: Use DimPlot with this named color vector
DimPlot(
  scRNA_epi_saline, 
  group.by = "CellType",   # Ensure "CellType" is the correct metadata for grouping
  reduction = "umap",
  label = F, 
  pt.size = 1, 
  cols = colors_use_larynx  # Named vector applies specific colors to specific identities
)

# Display the "alphabet" palette using show_palette() from ggsci
colors <- DiscretePalette_scCustomize(n = 24, palette = "alphabet")
# Use show_col to visualize the colors
show_col(colors)

#FOAOFFFF#0075DCFF#993F00FF#4C005CFF#191919FF
#005C31FF#2BCE48F#FFCC99F#808080FF#94FFB5FF
#8F7C00F#9DCC00FF#C20088FF#003380FF#FFA405FF
#FFA8BBF#426600FF#FF0010FF#5EF1F2FF#00998FFF
#E0FF66FF#740AFFFF#990000FF#FFFF80FF

Investigate the relationship between cluster identity and sample identity. Building a phylogenetic tree relating the ‘average’ cell from each group. This tree is estimated based on a distance matrix constructed in either gene expression space or PCA space.

#Per CELL TYPE
head(scRNA_epi_saline@meta.data)
##                                    orig.ident nCount_RNA nFeature_RNA
## AAACGCTCATTGGATC-1_1 all_compartments_saline1        706          485
## AAACGCTTCATACGAC-1_1 all_compartments_saline1      17139         4208
## AAACGCTTCGCGTGCA-1_1 all_compartments_saline1        677          495
## AAAGAACTCCCGAGGT-1_1 all_compartments_saline1       4073         1219
## AAAGAACTCCTTCAGC-1_1 all_compartments_saline1      42779         6137
## AAAGGATGTTGCACGC-1_1 all_compartments_saline1      11201         3608
##                      percent.mt       S.Score     G2M.Score Phase
## AAACGCTCATTGGATC-1_1   9.803922 -0.0392240323 -0.0610403871    G1
## AAACGCTTCATACGAC-1_1   6.981877 -0.0266484329 -0.0672562624    G1
## AAACGCTTCGCGTGCA-1_1   1.608187 -0.0397784353  0.0444191333   G2M
## AAAGAACTCCCGAGGT-1_1   6.619269  0.0207733279 -0.0007213987     S
## AAAGAACTCCTTCAGC-1_1   4.765234 -0.0006420594 -0.0539605231    G1
## AAAGGATGTTGCACGC-1_1   6.724705 -0.0372058950 -0.0355896740    G1
##                                     old.ident RNA_snn_res.0.5 seurat_clusters
## AAACGCTCATTGGATC-1_1 all_compartments_saline1               2              11
## AAACGCTTCATACGAC-1_1 all_compartments_saline1               2               2
## AAACGCTTCGCGTGCA-1_1 all_compartments_saline1               5              11
## AAAGAACTCCCGAGGT-1_1 all_compartments_saline1               3               4
## AAAGAACTCCTTCAGC-1_1 all_compartments_saline1               2               3
## AAAGGATGTTGCACGC-1_1 all_compartments_saline1               2               2
##                      pANN_0.1_0.01_83 DF.classifications_0.1_0.01_83 SampleName
## AAACGCTCATTGGATC-1_1          0.06250                        Singlet       <NA>
## AAACGCTTCATACGAC-1_1          0.09375                        Singlet       <NA>
## AAACGCTTCGCGTGCA-1_1          0.12500                        Singlet       <NA>
## AAAGAACTCCCGAGGT-1_1          0.03125                        Singlet       <NA>
## AAAGAACTCCTTCAGC-1_1          0.09375                        Singlet       <NA>
## AAAGGATGTTGCACGC-1_1          0.03125                        Singlet       <NA>
##                      Condition pANN_0.15_0.005_340
## AAACGCTCATTGGATC-1_1   Control                  NA
## AAACGCTTCATACGAC-1_1   Control                  NA
## AAACGCTTCGCGTGCA-1_1   Control                  NA
## AAAGAACTCCCGAGGT-1_1   Control                  NA
## AAAGAACTCCTTCAGC-1_1   Control                  NA
## AAAGGATGTTGCACGC-1_1   Control                  NA
##                      DF.classifications_0.15_0.005_340 integrated_snn_res.0.1
## AAACGCTCATTGGATC-1_1                              <NA>                      1
## AAACGCTTCATACGAC-1_1                              <NA>                      1
## AAACGCTTCGCGTGCA-1_1                              <NA>                      1
## AAAGAACTCCCGAGGT-1_1                              <NA>                      1
## AAAGAACTCCTTCAGC-1_1                              <NA>                      1
## AAAGGATGTTGCACGC-1_1                              <NA>                      1
##                      integrated_snn_res.0.8 integrated_snn_res.0.5
## AAACGCTCATTGGATC-1_1                     13                      1
## AAACGCTTCATACGAC-1_1                      2                      1
## AAACGCTTCGCGTGCA-1_1                     13                      1
## AAAGAACTCCCGAGGT-1_1                      3                      3
## AAAGAACTCCTTCAGC-1_1                      4                      4
## AAAGGATGTTGCACGC-1_1                      2                      1
##                      integrated_snn_res.1 integrated_snn_res.1.5
## AAACGCTCATTGGATC-1_1                   13                     13
## AAACGCTTCATACGAC-1_1                    2                      1
## AAACGCTTCGCGTGCA-1_1                   13                     13
## AAAGAACTCCCGAGGT-1_1                    3                     10
## AAAGAACTCCTTCAGC-1_1                    4                      2
## AAAGGATGTTGCACGC-1_1                    2                      1
##                      integrated_snn_res.1.2 integrated_snn_res.1.1
## AAACGCTCATTGGATC-1_1                     11                     13
## AAACGCTTCATACGAC-1_1                      2                      2
## AAACGCTTCGCGTGCA-1_1                     11                     13
## AAAGAACTCCCGAGGT-1_1                      4                      3
## AAAGAACTCCTTCAGC-1_1                      3                      4
## AAAGGATGTTGCACGC-1_1                      2                      2
##                           CellType           Region
## AAACGCTCATTGGATC-1_1 Club-proximal Tracheobronchial
## AAACGCTTCATACGAC-1_1 Club-proximal Tracheobronchial
## AAACGCTTCGCGTGCA-1_1 Club-proximal Tracheobronchial
## AAAGAACTCCCGAGGT-1_1   Club-distal Tracheobronchial
## AAAGAACTCCTTCAGC-1_1      Club-mid Tracheobronchial
## AAAGGATGTTGCACGC-1_1 Club-proximal Tracheobronchial
Idents(scRNA_epi_saline) <- scRNA_epi_saline$CellType
## extract meta data
membership <- scRNA_epi_saline@meta.data %>% as.data.table # the resulting membership object has one "row" per cell
## count the number of cells per unique combinations of "Condition" and "CellType"
membership[, .N, by = c("Condition", "CellType")]
##     Condition            CellType     N
##        <char>              <fctr> <int>
##  1:   Control       Club-proximal   532
##  2:   Control         Club-distal   335
##  3:   Control            Club-mid   336
##  4:   Control     Basal-v.trachea   527
##  5:   Control            Goblet-2   414
##  6:   Control     Basal-d.trachea   420
##  7:   Control            Goblet-1    88
##  8:   Control            Ciliated    48
##  9:   Control   Secretory-trachea   208
## 10:   Control    Parabasal-larynx   189
## 11:   Control         Serous-duct   148
## 12:   Control    Secretory-larynx    35
## 13:   Control        Serous-acini   116
## 14:   Control   Suprabasal-larynx    54
## 15:   Control        Basal-larynx   106
## 16:   Control Basal-myoepithelial    37
## 17:   Control                Tuft    19
## 18:   Control      Neuroendocrine     5
## with additional casting after the counting
membership[, .N, by = c("Condition", "CellType")] %>% dcast(., Condition ~ CellType, value.var = "N")
## Key: <Condition>
##    Condition Basal-v.trachea Basal-d.trachea Basal-larynx Parabasal-larynx
##       <char>           <int>           <int>        <int>            <int>
## 1:   Control             527             420          106              189
##    Basal-myoepithelial Suprabasal-larynx Secretory-trachea Secretory-larynx
##                  <int>             <int>             <int>            <int>
## 1:                  37                54               208               35
##    Club-proximal Club-mid Club-distal Goblet-1 Goblet-2 Serous-duct
##            <int>    <int>       <int>    <int>    <int>       <int>
## 1:           532      336         335       88      414         148
##    Serous-acini Neuroendocrine  Tuft Ciliated
##           <int>          <int> <int>    <int>
## 1:          116              5    19       48
dittoBarPlot(scRNA_epi_saline, "Condition", group.by = "CellType",
             main = "Cluster Composition",
             y.breaks = c(0,.1,.2,.3,.4,.5,.6,.7,.8,.9,1),
             xlab = NULL, # NULL = remove
             ylab = "% of total",
              color.panel = c("skyblue","brown"))

scRNA_epi_saline <- BuildClusterTree(scRNA_epi_saline, dims = 1:50)
PlotClusterTree(scRNA_epi_saline)

Total_epi_cell_count <- ggplot(scRNA_epi_saline@meta.data, aes(CellType)) +
  geom_bar(stat="count", fill='red', colour='black', size = 0.3, width = 0.8, key_glyph = draw_key_label, alpha=.6) +
labs(title = NULL, subtitle=NULL, x = NULL, y = 'Total Cell Count', caption = NULL, fill = NULL) + 
  theme(legend.text = element_text(size = 15, face=NULL)) +
  theme(plot.title=element_text(size=13, face="bold"), 
        axis.text.x=element_text(size=12, face="bold", angle = 40, hjust = 1), 
        axis.text.y=element_text(size=7),
        axis.title.x=element_text(size=8, face=NULL),
        axis.title.y=element_text(size =14)) +
        scale_y_continuous(expand = c(0, 0), breaks = seq(0, 1800, 100)) +
    coord_cartesian(ylim=c(0, 1800)) +
  theme(panel.grid = element_blank(), panel.background = element_blank(), axis.line = element_line(colour = 'black'), legend.position="right") 
VlnPlot_scCustom(scRNA_epi_saline, features = "nCount_RNA", split.by = "Condition", group.by = "CellType", pt.size = 0)

#Explore DEG with advanced visualization

#Define assay
DefaultAssay(scRNA_epi_saline) <- "RNA"
#SMG subpopulations
name <- "SMG_populations"
FeaturePlot_scCustom(scRNA_epi_saline, 
            features = c("Dmbt1","Acta2","Nkx3-1","Muc5b","Tff2","Lipf","Slc34a2","Bpifb5"), 
            order=TRUE, pt.size=1, colors_use= c("#E4E1E3FF","#F6222EFF"), na_cutoff =.5, num_columns = 4)

FeaturePlot_scCustom(scRNA_epi_saline, 
            features = c("Dmbt1"), 
            order=TRUE, pt.size=2, colors_use= c("#E4E1E3FF","#F6222EFF"), na_cutoff =2)

FeaturePlot_scCustom(scRNA_epi_saline, 
            features = c("Acta2"), 
            order=TRUE, pt.size=2, colors_use= c("#E4E1E3FF","#F6222EFF"), na_cutoff =1)

FeaturePlot_scCustom(scRNA_epi_saline, 
            features = c("Nkx3-1"), 
            order=TRUE, pt.size=2, colors_use= c("#E4E1E3FF","#F6222EFF"), na_cutoff =.5)

FeaturePlot_scCustom(scRNA_epi_saline, 
            features = c("Muc5b"), 
            order=TRUE, pt.size=2, colors_use= c("#E4E1E3FF","#F6222EFF"), na_cutoff =3)

FeaturePlot_scCustom(scRNA_epi_saline, 
            features = c("Tff2"), 
            order=TRUE, pt.size=2, colors_use= c("#E4E1E3FF","#F6222EFF"), na_cutoff =3)

FeaturePlot_scCustom(scRNA_epi_saline, 
            features = c("Lipf"), 
            order=TRUE, pt.size=2, colors_use= c("#E4E1E3FF","#F6222EFF"), na_cutoff =4)

FeaturePlot_scCustom(scRNA_epi_saline, 
            features = c("Slc34a2"), 
            order=TRUE, pt.size=2, colors_use= c("#E4E1E3FF","#F6222EFF"), na_cutoff =.5)

FeaturePlot_scCustom(scRNA_epi_saline, 
            features = c("Bpifb5"), 
            order=TRUE, pt.size=2, colors_use= c("#E4E1E3FF","#F6222EFF"), na_cutoff =.5)

FeaturePlot_scCustom(scRNA_epi_saline, 
            features = c("Nkx3-1","Sox9","Bpifb2","Azgp1"), 
            order=TRUE, pt.size=2, colors_use= c("#E4E1E3FF","#F6222EFF"), na_cutoff =.5)

#Basal-to-luminal markers
FeaturePlot_scCustom(scRNA_epi_saline, 
            features = c("Igfbp2"), 
            order=TRUE,pt.size=2, colors_use= c("#E4E1E3FF","#F6222EFF"), na_cutoff =1.5)

FeaturePlot_scCustom(scRNA_epi_saline, 
            features = c("Tmprss11a"), 
            order=TRUE,pt.size=2, colors_use= c("#E4E1E3FF","#F6222EFF"), na_cutoff =1.5)

FeaturePlot_scCustom(scRNA_epi_saline, 
            features = c("Krt13"), 
            order=TRUE,pt.size=2, colors_use= c("#E4E1E3FF","#F6222EFF"), na_cutoff =4)

#Regional tissue markers
FeaturePlot_scCustom(scRNA_epi_saline, 
            features = c("Tmprss11a"), 
            order=TRUE,pt.size=2, colors_use= c("#E4E1E3FF","#F6222EFF"), na_cutoff =1.5)

FeaturePlot_scCustom(scRNA_epi_saline, 
            features = c("Nkx2-1"), 
            order=TRUE,pt.size=2, colors_use= c("#E4E1E3FF","#F6222EFF"), na_cutoff =.5)

FeaturePlot_scCustom(scRNA_epi_saline, 
            features = c("Dmbt1"), 
            order=TRUE,pt.size=2, colors_use= c("#E4E1E3FF","#F6222EFF"), na_cutoff =2)

FeaturePlot_scCustom(scRNA_epi_saline, 
            features = c("Sox9"), 
            order=TRUE,pt.size=2, colors_use= c("#E4E1E3FF","#F6222EFF"), na_cutoff =.5)

#Club markers
FeaturePlot_scCustom(scRNA_epi_saline, 
            features = c("Scgb1a1"), 
            order=TRUE,pt.size=1, colors_use= c("#E4E1E3FF","#F6222EFF"), na_cutoff =4)

FeaturePlot_scCustom(scRNA_epi_saline, 
            features = c("Scgb3a2"), 
            order=TRUE,pt.size=1, colors_use= c("#E4E1E3FF","#F6222EFF"), na_cutoff =4)

#Keratin markers
FeaturePlot_scCustom(scRNA_epi_saline, 
            features = c('Krt1', 'Krt2', 'Krt3', 'Krt4', 'Krt5', 'Krt6a', 'Krt6b', 'Krt6c', 'Krt7', 'Krt8', 'Krt9', 'Krt10', 'Krt12', 'Krt13', 'Krt14', 'Krt15', 'Krt16', 'Krt17', 'Krt18', 'Krt19', 'Krt20', 'Krt21', 'Krt23', 'Krt24', 'Krt25', 'Krt26', 'Krt27', 'Krt28', 'Krt31', 'Krt32', 'Krt33a', 'Krt33b', 'Krt34', 'Krt35', 'Krt36', 'Krt37', 'Krt38', 'Krt39', 'Krt40', 'Krt71', 'Krt72', 'Krt73', 'Krt74', 'Krt75', 'Krt76', 'Krt77', 'Krt78', 'Krt79', 'Krt80', 'Krt81', 'Krt82', 'Krt83', 'Krt84', 'Krt85', 'Krt86', 'Krt87', 'Krt88', 'Krt89', 'Krt90', 'Krt93', 'Krt94', 'Krt95', 'Krt96'), 
            colors_use= c("#E4E1E3FF","#F6222EFF"), na_cutoff =.5, num_columns = 3)

DotPlot_scCustom(scRNA_epi_saline, 
            features = c('Krt1', 'Krt2', 'Krt3', 'Krt4', 'Krt5', 'Krt6a', 'Krt6b', 'Krt6c', 'Krt7', 'Krt8', 'Krt9', 'Krt10', 'Krt12', 'Krt13', 'Krt14', 'Krt15', 'Krt16', 'Krt17', 'Krt18', 'Krt19', 'Krt20', 'Krt21', 'Krt23', 'Krt24', 'Krt25', 'Krt26', 'Krt27', 'Krt28', 'Krt31', 'Krt32', 'Krt33a', 'Krt33b', 'Krt34', 'Krt35', 'Krt36', 'Krt37', 'Krt38', 'Krt39', 'Krt40', 'Krt71', 'Krt72', 'Krt73', 'Krt74', 'Krt75', 'Krt76', 'Krt77', 'Krt78', 'Krt79', 'Krt80', 'Krt81', 'Krt82', 'Krt83', 'Krt84', 'Krt85', 'Krt86', 'Krt87', 'Krt88', 'Krt89', 'Krt90', 'Krt93', 'Krt94', 'Krt95', 'Krt96'), x_lab_rotate = TRUE,
            colors_use= c("#E4E1E3FF","#F6222EFF")) #63 known KRTs

Session information

sessionInfo()
## R version 4.3.3 (2024-02-29)
## Platform: aarch64-apple-darwin20 (64-bit)
## Running under: macOS 15.2
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/lib/libRblas.0.dylib 
## LAPACK: /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/lib/libRlapack.dylib;  LAPACK version 3.11.0
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## time zone: America/Los_Angeles
## tzcode source: internal
## 
## attached base packages:
## [1] stats4    grid      stats     graphics  grDevices utils     datasets 
## [8] methods   base     
## 
## other attached packages:
##  [1] EnhancedVolcano_1.20.0 ggrepel_0.9.6          biomaRt_2.58.2        
##  [4] org.Mm.eg.db_3.18.0    enrichR_3.2            topGO_2.54.0          
##  [7] SparseM_1.84-2         GO.db_3.18.0           AnnotationDbi_1.64.1  
## [10] IRanges_2.36.0         S4Vectors_0.40.2       Biobase_2.62.0        
## [13] graph_1.80.0           BiocGenerics_0.48.1    igraph_2.1.2          
## [16] paletteer_1.6.0        qs_0.27.2              scCustomize_3.0.1     
## [19] data.table_1.16.4      dittoSeq_1.14.3        viridis_0.6.5         
## [22] viridisLite_0.4.2      magrittr_2.0.3         ComplexHeatmap_2.18.0 
## [25] HGNChelper_0.8.15      scales_1.3.0           ggsci_3.2.0           
## [28] ggplot2_3.5.1          dplyr_1.1.4            tidyr_1.3.1           
## [31] kableExtra_1.4.0       Seurat_5.1.0           SeuratObject_5.0.2    
## [34] sp_2.1-4              
## 
## loaded via a namespace (and not attached):
##   [1] matrixStats_1.5.0           spatstat.sparse_3.1-0      
##   [3] bitops_1.0-9                lubridate_1.9.4            
##   [5] httr_1.4.7                  RColorBrewer_1.1-3         
##   [7] doParallel_1.0.17           tools_4.3.3                
##   [9] sctransform_0.4.1           R6_2.5.1                   
##  [11] lazyeval_0.2.2              uwot_0.1.16                
##  [13] GetoptLong_1.0.5            withr_3.0.2                
##  [15] prettyunits_1.2.0           gridExtra_2.3              
##  [17] progressr_0.15.1            cli_3.6.3                  
##  [19] Cairo_1.6-2                 spatstat.explore_3.3-4     
##  [21] fastDummies_1.7.4           prismatic_1.1.2            
##  [23] labeling_0.4.3              sass_0.4.9                 
##  [25] spatstat.data_3.1-4         ggridges_0.5.6             
##  [27] pbapply_1.7-2               systemfonts_1.1.0          
##  [29] svglite_2.1.3               parallelly_1.41.0          
##  [31] WriteXLS_6.7.0              limma_3.58.1               
##  [33] rstudioapi_0.17.1           RSQLite_2.3.9              
##  [35] generics_0.1.3              shape_1.4.6.1              
##  [37] RApiSerialize_0.1.4         ica_1.0-3                  
##  [39] spatstat.random_3.3-2       Matrix_1.6-5               
##  [41] ggbeeswarm_0.7.2            abind_1.4-8                
##  [43] lifecycle_1.0.4             yaml_2.3.10                
##  [45] snakecase_0.11.1            SummarizedExperiment_1.32.0
##  [47] BiocFileCache_2.10.2        SparseArray_1.2.4          
##  [49] Rtsne_0.17                  blob_1.2.4                 
##  [51] promises_1.3.2              crayon_1.5.3               
##  [53] miniUI_0.1.1.1              lattice_0.22-6             
##  [55] cowplot_1.1.3               KEGGREST_1.42.0            
##  [57] pillar_1.10.1               knitr_1.49                 
##  [59] GenomicRanges_1.54.1        rjson_0.2.23               
##  [61] future.apply_1.11.3         codetools_0.2-20           
##  [63] leiden_0.4.3.1              glue_1.8.0                 
##  [65] spatstat.univar_3.1-1       vctrs_0.6.5                
##  [67] png_0.1-8                   spam_2.11-0                
##  [69] gtable_0.3.6                rematch2_2.1.2             
##  [71] cachem_1.1.0                xfun_0.50                  
##  [73] S4Arrays_1.2.1              mime_0.12                  
##  [75] survival_3.8-3              SingleCellExperiment_1.24.0
##  [77] pheatmap_1.0.12             iterators_1.0.14           
##  [79] statmod_1.5.0               fitdistrplus_1.2-2         
##  [81] ROCR_1.0-11                 nlme_3.1-166               
##  [83] bit64_4.6.0-1               filelock_1.0.3             
##  [85] progress_1.2.3              RcppAnnoy_0.0.22           
##  [87] GenomeInfoDb_1.38.8         bslib_0.8.0                
##  [89] irlba_2.3.5.1               vipor_0.4.7                
##  [91] KernSmooth_2.23-26          splitstackshape_1.4.8      
##  [93] colorspace_2.1-1            DBI_1.2.3                  
##  [95] ggrastr_1.0.2               tidyselect_1.2.1           
##  [97] curl_6.1.0                  bit_4.5.0.1                
##  [99] compiler_4.3.3              xml2_1.3.6                 
## [101] DelayedArray_0.28.0         plotly_4.10.4              
## [103] stringfish_0.16.0           lmtest_0.9-40              
## [105] rappdirs_0.3.3              stringr_1.5.1              
## [107] digest_0.6.37               goftest_1.2-3              
## [109] presto_1.0.0                spatstat.utils_3.1-2       
## [111] rmarkdown_2.29              XVector_0.42.0             
## [113] htmltools_0.5.8.1           pkgconfig_2.0.3            
## [115] MatrixGenerics_1.14.0       dbplyr_2.5.0               
## [117] fastmap_1.2.0               rlang_1.1.4                
## [119] GlobalOptions_0.1.2         htmlwidgets_1.6.4          
## [121] shiny_1.10.0                farver_2.1.2               
## [123] jquerylib_0.1.4             zoo_1.8-12                 
## [125] jsonlite_1.8.9              RCurl_1.98-1.16            
## [127] GenomeInfoDbData_1.2.11     dotCall64_1.2              
## [129] patchwork_1.3.0             munsell_0.5.1              
## [131] Rcpp_1.0.14                 ape_5.8-1                  
## [133] reticulate_1.40.0           stringi_1.8.4              
## [135] zlibbioc_1.48.2             MASS_7.3-60.0.1            
## [137] plyr_1.8.9                  parallel_4.3.3             
## [139] listenv_0.9.1               forcats_1.0.0              
## [141] deldir_2.0-4                Biostrings_2.70.3          
## [143] splines_4.3.3               tensor_1.5                 
## [145] hms_1.1.3                   circlize_0.4.16            
## [147] spatstat.geom_3.3-4         RcppHNSW_0.6.0             
## [149] reshape2_1.4.4              XML_3.99-0.18              
## [151] evaluate_1.0.3              RcppParallel_5.1.9         
## [153] ggprism_1.0.5               foreach_1.5.2              
## [155] httpuv_1.6.15               RANN_2.6.2                 
## [157] purrr_1.0.2                 polyclip_1.10-7            
## [159] future_1.34.0               clue_0.3-66                
## [161] scattermore_1.2             janitor_2.2.1              
## [163] xtable_1.8-4                RSpectra_0.16-2            
## [165] later_1.4.1                 tibble_3.2.1               
## [167] memoise_2.0.1               beeswarm_0.4.0             
## [169] cluster_2.1.8               timechange_0.3.0           
## [171] globals_0.16.3